
class Student : public Distribution 
{
        public:
        Student(){}
        Student(int n,int nu, mat m, mat s)
        {
                _n=n;
                _m=m;
                _nu=nu;
                _s=s;//chol here
                _cs=chol(s);
                _is=inv(s);
                boost::random::student_t_distribution<> N(nu);
                //dn= new boost::math::student_t_distribution<>(nu); 
               foo1=tgamma((double)(_nu+_n)/2);
               bar1=tgamma((double)(_nu)/2);
		td=new boost::math::students_t_distribution<>(nu);
	       	rt= new RandomG::Random<boost::random::student_t_distribution<> >(N);
		boost::random::uniform_01<> U;
		runif=new RandomG::Random<boost::random::uniform_01<> >(U);
        }
        virtual ~Student()
        {
                delete rt;
                delete td;
               // delete dn;
        }
	
	mat scrambled(int M){
		mat X(M,_n);
		double *N=new double[M*_n];
		sobol::Owenscramble O;
		O(M,_n,0,N);
		for(int i=0;i<(M*_n);i++)
		{
			if(N[i]<0.0000000001)
			{
				N[i]=(*runif)();// TEmporary problem with first value 4 vales of the sequence
			}
			N[i]=quantile(*td,N[i]);

		}
		X=vect2mat<double>(N,M,_n);
		mat temp;
		for( int i=0;i<M;i++)
		{
			temp=(_m+(_cs).t()*X(i,span::all).t());
			//X(i,span::all)=(this->r(1)).t();//temp.t();
		}

		return X;
	}
        mat r(int j)
        {
                mat X(_n,1);
                for(int i=0;i<_n;i++)
                {
                        X(i,0)=(*rt)();
                }
                X=_m+(_cs).t()*X;
                return X;
        }


        double  d(mat x, int blog=0,int kk=0){
        
                double res=0;//
		if(blog)
		{
				mat foo=x-_m;
				mat bar=1+((double)1/_nu)*foo.t()*_is*foo;
				//cout << "//" << bar;
				res=log(foo1)-log(bar1)-0.5*_n*log(_nu*PI);
				res+=-0.5*log(det(_s))-0.5*(_nu+_n)*log(as_scalar(bar));
				if(res<_mINF)
				{
					cout << " problem !!!!";
					res=_mINF;
				}	
		}else{
			mat foo=x-_m;
			mat bar=1+((double)1/_nu)*foo.t()*_is*foo;
			cout << bar1;

			double temp=pow(_nu*PI,_n)*det(_s)*pow(bar(0,0),(_nu+_n));
			cout << "temp " <<temp;
			res=foo1/(bar1*sqrt(temp));
			if(res==0)
			{
				cout << " problem !!!!";
				res=0.000001;
			}
		}
		return res;
        }
        mat Var(mat X){return cov(X);}//((_nu-2)/_nu)*cov(X);}
	void Set_n(int n){     _n=n;   }
        void Set_m(mat m){      _m=m;   }
        mat Get_m(void){      return _m;   }
        mat Get_s(void){      return _s;   }
        void Set_s(mat s){  _cs=chol(s);_is=inv(s);_s=s;  }
        private:
        int _n;
        mat _m;
	double foo1;
	double bar1;
        mat _s;
        mat _is;
        mat _cs;
	int _nu;
        RandomG::Random<boost::random::student_t_distribution<> > *rt;
	RandomG::Random<boost::random::uniform_01<> > *runif;
 	boost::math::students_t_distribution<> *td;
	//boost::math::normal_distribution<> *dn;

};
